ANALYSIS ON GLOBAL WARMING AND THE PRODUCTION OF C02
IMPORT LIBRARIES
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import plotly
import plotly.graph_objs as go
import plotly.tools as tls
import seaborn as sns
import time
import warnings
import cartopy.crs as ccrs
warnings.filterwarnings('ignore')
import os
import plotly.express as px
import scipy
import datetime
from datetime import datetime
import pickle
import matplotlib.dates as mdates
import xarray as xr
FUNCTIONS
def fetch_year(date):
return date.split('-')[0]
def get_season (month):
if month>=3 and month<=5:
return 'spring'
elif month>=6 and month<=8:
return 'summer'
elif month>=9 and month<=11:
return 'autumn'
else:
return 'winter'
IMPORT DATASETS
The dataasets are available at following links: -https://www.kaggle.com/datasets/berkeleyearth/climate-change-earth-surface-temperature-data -https://github.com/owid/co2-data -The models are some of the CMIP6 models that we have seen during the lessons
#os.chdir('C:\\Users\\User\\Documents\\Magistrale_SSE\\lab_Geographic info systems\\_PROGETTOFINALE2')
#os.getcwd()
global_temp = pd.read_csv('GlobalTemperatures.csv')
global_temp_country=pd.read_csv('GlobalLandTemperaturesByCountry.csv')
co2=pd.read_csv('owid-co2-data.csv')
os.getcwd()
'C:\\Users\\User\\Documents\\Magistrale_SSE\\lab_Geographic info systems\\_PROGETTOFINALE2'
modfile1='ts_Amon_IPSL-CM6A-LR_historical_r1i1p1f1_gr_185001-201412.nc'
modfile2='ts_Amon_MIROC6_historical_r1i1p1f1_gn_185001-201412.nc'
modfile3='ts_Amon_MRI-ESM2-0_historical_r1i1p1f1_gn_185001-201412.nc'
AT THE BEGININNING WE ANALYZE THE AVERAGE LAND TEMPERATURE IN THE COURSE OF THE YEARS, THEN WE ANALYZE THE PRODUCTION OF CO2
global_temp['years']=global_temp['dt'].apply(fetch_year)
global_temp=global_temp[global_temp['years']>='1850']
global_temp.groupby('years').agg({'LandAverageTemperature':'mean','LandAverageTemperatureUncertainty':'mean'})
| LandAverageTemperature | LandAverageTemperatureUncertainty | |
|---|---|---|
| years | ||
| 1850 | 7.900667 | 0.876417 |
| 1851 | 8.178583 | 0.881917 |
| 1852 | 8.100167 | 0.918250 |
| 1853 | 8.041833 | 0.835000 |
| 1854 | 8.210500 | 0.825667 |
| ... | ... | ... |
| 2011 | 9.516000 | 0.082000 |
| 2012 | 9.507333 | 0.083417 |
| 2013 | 9.606500 | 0.097667 |
| 2014 | 9.570667 | 0.090167 |
| 2015 | 9.831000 | 0.092167 |
166 rows × 2 columns
data=global_temp.groupby('years').agg({'LandAverageTemperature':'mean','LandAverageTemperatureUncertainty':'mean'}).reset_index()
#creiamo l'estremo superiore e l'estremo inferiore del nostro intervallo
data['Uncertainty top']=data['LandAverageTemperature']+data['LandAverageTemperatureUncertainty']
data['Uncertainty bottom']=data['LandAverageTemperature']-data['LandAverageTemperatureUncertainty']
fig=px.line(data,x='years',y=['LandAverageTemperature',
'Uncertainty top', 'Uncertainty bottom'],title='Average Land Tmeperature in World')
fig.show()
Let's focus ourselves on the seasons
global_temp['dt']=pd.to_datetime(global_temp['dt'])
global_temp['month']=global_temp['dt'].dt.month
global_temp['season']=global_temp['month'].apply(get_season)
years=global_temp['years'].unique()
spring_temps=[]
summer_temps=[]
autumn_temps=[]
winter_temps=[]
for year in years:
current_year=global_temp[global_temp['years']==year]
spring_temps.append(current_year[current_year['season']=='spring']['LandAverageTemperature'].mean())
summer_temps.append(current_year[current_year['season']=='summer']['LandAverageTemperature'].mean())
autumn_temps.append(current_year[current_year['season']=='autumn']['LandAverageTemperature'].mean())
winter_temps.append(current_year[current_year['season']=='winter']['LandAverageTemperature'].mean())
#Now let's make a dataframe on season
season=pd.DataFrame()
season['year']=years
season['spring_temp']=spring_temps
season['summer_temp']=summer_temps
season['autumn_temp']=autumn_temps
season['winter_temp']=winter_temps
fig=px.line(season,x='year',y=['spring_temp', 'summer_temp', 'autumn_temp', 'winter_temp'],title='Average temperature in Each season')
fig.show()
From the charts we can see, that there is global warming nowadays. The average temperature of Earth surface has the highest value in the last three centuries. The fastest temperature growth occurred in the last 30 years. This charts also have confidence intervals, which shows that measurement of temperature has become more accurate in the last few years.
FOCUS ON CO2 EMISSION
co2_inter=co2[['country','year','co2','gdp']]
co2_inter=co2_inter[co2_inter['country']!='World']
co2_agg=co2_inter.groupby('year').agg({'co2':'mean'}).reset_index()
fig=px.line(co2_agg,x='year',y='co2',title='Mean of co2 during the years')
fig.show()
from this graphic we can see that has been an augmentation of CO2 in the atmosphere from the 1950, let's do a test between the mean of co2 and the land average temperature
global_agg=global_temp.groupby('years').agg({'LandAverageTemperature':'mean'}).reset_index() #temperatura media per anno
global_agg=pd.DataFrame(global_agg)
co2_agg=co2_agg[co2_agg['year']<=2015]
co2_agg=co2_agg[co2_agg['year']>=1850]
from scipy import stats
scipy.stats.pearsonr(co2_agg['co2'], global_agg['LandAverageTemperature'])
PearsonRResult(statistic=0.8751956131048915, pvalue=1.3540117283304936e-53)
We refuse the hypothesis that the variables are not correlated,so CO2 is correlated with the Land Average Temperature
Let's now see if gdp is correlated with the augmentation of CO2
co2_gdp=co2_inter.groupby('year').agg({'gdp':'sum'}).reset_index()
co2_agg=co2_agg[co2_agg['year']<=2015]
co2_agg=co2_agg[co2_agg['year']>=1850]
co2_gdp=co2_gdp[co2_gdp['year']<=2015]
co2_gdp=co2_gdp[co2_gdp['year']>=1850]
from scipy import stats
scipy.stats.pearsonr(co2_agg['co2'], co2_gdp['gdp'])
PearsonRResult(statistic=0.9395288360896425, pvalue=3.152489761450195e-78)
In statistics correlation is not synonymous of causality, but from this analysis we can see that there are bonds beteween the global average land co2 and the gdp of a country
Let's make a comparison between model ensemble from CMIP6 (this models simulate the augmentation of CO2 in the course of the years) and the current situation of the Average Land Temperature
d1 = xr.open_dataset(modfile1)
d2 = xr.open_dataset(modfile2)
d3 = xr.open_dataset(modfile3)
d1
<xarray.Dataset>
Dimensions: (lat: 143, lon: 144, time: 1980, axis_nbounds: 2)
Coordinates:
* lat (lat) float32 -90.0 -88.73 -87.46 -86.2 ... 87.46 88.73 90.0
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* time (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:...
Dimensions without coordinates: axis_nbounds
Data variables:
time_bounds (time, axis_nbounds) datetime64[ns] 1850-01-01 ... 2015-01-01
ts (time, lat, lon) float32 ...
Attributes: (12/51)
Conventions: CF-1.7 CMIP-6.2
creation_date: 2018-07-11T07:36:33Z
tracking_id: hdl:21.14100/3168f5b1-bf0a-4aec-931f-73c9d0034a45
description: CMIP6 historical
title: IPSL-CM6A-LR model output prepared for CMIP6 / CM...
activity_id: CMIP
... ...
name: /ccc/work/cont003/gencmip6/p86caub/IGCM_OUT/IPSLC...
further_info_url: https://furtherinfo.es-doc.org/CMIP6.IPSL.IPSL-CM...
variant_label: r1i1p1f1
realization_index: 1
history: Sat Dec 1 12:17:27 2018: ncatted -O -a realizati...
NCO: "4.6.0"d2
<xarray.Dataset>
Dimensions: (time: 1980, bnds: 2, lat: 128, lon: 256)
Coordinates:
* time (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00
* lat (lat) float64 -88.93 -87.54 -86.14 -84.74 ... 86.14 87.54 88.93
* lon (lon) float64 0.0 1.406 2.812 4.219 ... 354.4 355.8 357.2 358.6
Dimensions without coordinates: bnds
Data variables:
time_bnds (time, bnds) datetime64[ns] 1850-01-01 1850-02-01 ... 2015-01-01
lat_bnds (lat, bnds) float64 -90.0 -88.28 -88.28 ... 88.28 88.28 90.0
lon_bnds (lon, bnds) float64 -0.7031 0.7031 0.7031 ... 357.9 357.9 359.3
ts (time, lat, lon) float32 ...
Attributes: (12/45)
Conventions: CF-1.7 CMIP-6.2
activity_id: CMIP
branch_method: standard
branch_time_in_child: 0.0
branch_time_in_parent: 0.0
creation_date: 2018-11-30T16:15:09Z
... ...
variable_id: ts
variant_label: r1i1p1f1
license: CMIP6 model data produced by MIROC is licensed un...
cmor_version: 3.3.2
tracking_id: hdl:21.14100/24645cf7-2812-40bc-a320-cfc906678afe
NCO: netCDF Operators version 4.7.6 (Homepage = http:/...d3
<xarray.Dataset>
Dimensions: (time: 1980, bnds: 2, lat: 160, lon: 320)
Coordinates:
* time (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00
* lat (lat) float64 -89.14 -88.03 -86.91 -85.79 ... 86.91 88.03 89.14
* lon (lon) float64 0.0 1.125 2.25 3.375 ... 355.5 356.6 357.8 358.9
Dimensions without coordinates: bnds
Data variables:
time_bnds (time, bnds) datetime64[ns] 1850-01-01 1850-02-01 ... 2015-01-01
lat_bnds (lat, bnds) float64 -90.0 -88.59 -88.59 ... 88.59 88.59 90.0
lon_bnds (lon, bnds) float64 -0.5625 0.5625 0.5625 ... 358.3 358.3 359.4
ts (time, lat, lon) float32 ...
Attributes: (12/44)
Conventions: CF-1.7 CMIP-6.2
activity_id: CMIP
branch_method: standard
branch_time_in_child: 0.0
branch_time_in_parent: 0.0
creation_date: 2019-02-20T02:26:52Z
... ...
title: MRI-ESM2-0 output prepared for CMIP6
variable_id: ts
variant_label: r1i1p1f1
license: CMIP6 model data produced by MRI is licensed unde...
cmor_version: 3.4.0
tracking_id: hdl:21.14100/83b81c22-3ba4-45c0-9e9b-3f94038f7cb7#let's visualize some models, visualizzazione dell'anno 2014
p0 = d1.ts.isel(time=-1).plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet')
p0.axes.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x2b39b6dbdf0>
p1=d2.ts.isel(time=-1).plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet')
p1.axes.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x2b39b923340>
p2=d3.ts.isel(time=-1).plot(transform=ccrs.PlateCarree(),subplot_kws={'projection': ccrs.Robinson()},
cmap='jet')
p2.axes.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x2b39b7e1ff0>
#calculate normals
n1=d1.ts.sel(time=slice("1961-01-16T12:00:00.000000000","1990-12-16T12:00:00.000000000")
).groupby("time.month").mean()
n2=d2.ts.sel(time=slice("1961-01-16T12:00:00.000000000","1990-12-16T12:00:00.000000000")
).groupby("time.month").mean()
n3=d3.ts.sel(time=slice("1961-01-16T12:00:00.000000000","1990-12-16T12:00:00.000000000")
).groupby("time.month").mean()
# Second,calculate temperature anomalies
anom1 = d1.ts.groupby("time.month")-n1
anom2 = d2.ts.groupby("time.month")-n2
anom3 = d3.ts.groupby("time.month")-n3
# Weights
weights1 = np.cos(np.deg2rad(d1.lat))
weights2 = np.cos(np.deg2rad(d2.lat))
weights3 = np.cos(np.deg2rad(d3.lat))
#Aggregation
enst = xr.Dataset()
enst['ts1'] = anom1.weighted(weights1).mean(("lon", "lat"),keep_attrs=True).resample(time="Y").mean()
enst['ts2'] = anom2.weighted(weights2).mean(("lon", "lat")).resample(time="Y").mean()
enst['ts3'] = anom3.weighted(weights3).mean(("lon", "lat")).resample(time="Y").mean()
Let's make a comparison between the real situation and the situation predicted from the model
plt.plot(enst.ts1)
plt.plot(enst.ts2)
plt.plot(enst.ts3)
[<matplotlib.lines.Line2D at 0x2b395ba4d60>]
# Fourth, we can calculate the ensemble mean time series of global mean temperature anomalies
# Calculate the standard deviation of the 3 models
ens_mean = enst.to_array(dim='new').mean('new')
ens_std = enst.to_array(dim='new2').std('new2')
enst['ens_mean'] = ens_mean
enst['ens_std'] = ens_std
data=global_temp.groupby('years').agg({'LandAverageTemperature':'mean','LandAverageTemperatureUncertainty':'mean'}).reset_index()
obs=data[data['years']>='1850']
obs=obs[obs['years']<'2015']
ref=obs[obs['years']>='1961']
ref=ref[ref['years']<='1990']
Let's calculate anomalies for the temperatures, from our dataset global_temp_country and confront them with the anomalies of the model ensemble
media=ref.mean()['LandAverageTemperature']
#CALCOLO ANOMALIE
obs['temp2']=obs['LandAverageTemperature']-media
t=enst.ens_mean
p=enst.time
t=pd.DataFrame(t,p)
t=t.reset_index()
t['data']=pd.to_datetime(t['index'],format='%Y%m%d')
t['year']=pd.DatetimeIndex(t['data']).year
t.rename(columns={0:'temp'},inplace=True)
t
| index | temp | data | year | |
|---|---|---|---|---|
| 0 | 1850-12-31 | -0.259619 | 1850-12-31 | 1850 |
| 1 | 1851-12-31 | -0.216288 | 1851-12-31 | 1851 |
| 2 | 1852-12-31 | -0.240269 | 1852-12-31 | 1852 |
| 3 | 1853-12-31 | -0.242512 | 1853-12-31 | 1853 |
| 4 | 1854-12-31 | -0.147570 | 1854-12-31 | 1854 |
| ... | ... | ... | ... | ... |
| 160 | 2010-12-31 | 0.526859 | 2010-12-31 | 2010 |
| 161 | 2011-12-31 | 0.410785 | 2011-12-31 | 2011 |
| 162 | 2012-12-31 | 0.547878 | 2012-12-31 | 2012 |
| 163 | 2013-12-31 | 0.628964 | 2013-12-31 | 2013 |
| 164 | 2014-12-31 | 0.558643 | 2014-12-31 | 2014 |
165 rows × 4 columns
obs=obs.reset_index()
t['temp2']=obs['temp2']
t.rename(columns={'temp2':'reality','temp':'simulated'},inplace=True)
fig=px.line(t,x='year',y=['simulated','reality'],title='Anomalies in Temperature in reality vs simulated')
fig.show()
From this graphic we can see that the Global Land Temperature is increased more than the models have predicted, it seems interesting to confront the growing trend from the 1950. In this year, from the graphic of CO2 in the course of years that we have shown before, we can see that there has been an increase of CO2 that coincides with the growth of trend in this graphic. Our ensemble of models has predcited quite well the average temperature